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Understanding ballistic phonon transport effects in transient thermorefiectance experiments and 
explaining the observed deviations from classical theory remains a challenge. Diffusion equations are 
simple and computationally efficient but are widely believed to break down when the characteristic 
length scale is similar or less than the phonon mean-free-path. Building on our prior work, we 
demonstrate how well-known diffusion equations, namely the hyperbolic heat equation and the 
Cattaneo equation, can be used to model ballistic phonon effects in frequency-dependent periodic 
steady-state thermal transport. Our analytical solutions are found to compare excellently to rigorous 
numerical results of the phonon Boltzmann transport equation. The correct physical boundary 
conditions can be different from those traditionally used and are paramount for accurately capturing 
ballistic effects. To illustrate the technique, we consider a simple model problem using two different, 
commonly-used heating conditions. We demonstrate how this framework can easily handle detailed 
material properties, by considering the case of bulk silicon using a full phonon dispersion and 
mean-free-path distribution. This physically transparent approach provides clear insights into the 
nonequilibrium physics of quasi-ballistic phonon transport and its impact on thermal transport 
properties. 


I. INTRODUCTION 

Recent transient experiments probing the thermal 
transport properties of materials on short length- and/or 
time-scales have reported deviations from expected clas¬ 
sical theory, which often corresponds to a reduction in 
extracted thermal conductivity This has largely 

been attributed to quasi-ballistic phonon transport, aris¬ 
ing when the phonon mean-free-path (MFP) is similar to 
or greater than the characteristic length scale in the ex¬ 
periment. Traditional diffusion heat equations are com¬ 
monly relied upon for analyzing raw data and extracting 
thermal properties, but the classical heat equations are 
widely believed to break down under conditions of non- 
diffusive transport. This paper addresses the need for 
fast and accurate techniques to analyze transient ther¬ 
mal measurements. 

To capture ballistic behavior, theoretical efforts have 
largely focused on using rigorous approaches such as the 
phonon Boltzmann transport equation (BTE). Such de¬ 
tailed numerical studies (THI^ can be too computation¬ 
ally demanding for the routine analysis of experiments. 
The phonon BTE can, however, also provide a starting 
point for deriving simple models to elucidate the rela¬ 
tionship between the intrinsic phonon properties of a 
material (e.g. phonon MFP’) and the measured thermal 
transport characteristics [l3l - l^ . These simple models 
are often problem-specific - assuming simplified geome¬ 
tries and material structures. Moreover they are not 


‘Electronic address: jmaassen@dal.ca 


straightfowardly compatible with traditional analysis ap¬ 
proaches, and must often be used as post-processing tools 
to analyze the extracted thermal properties. Ideally one 
would like to have an approach that captures ballistic 
effects, but that can also be applied to a wide class of 
problems, readily handle material structures similar to 
the experimental setup and that can be used to analyze 
raw data. Such a technique is described in this paper. 

In our previous work, we showed that steady-state and 
transient diffusion equations can capture ballistic phonon 
effects as long as the correct physical boundary condi¬ 
tions are used. When properly implemented, diffusion 
equations provide good agreement with the phonon BTE 
1^ . In this paper, we extend previous work to the 
periodic steady-state case and analyze ballistic effects in 
model transient thermorefiectance experiments using dif¬ 
fusion equations. We compare diffusion equation solu¬ 
tions to recently-reported rigorous numerical results of 
the phonon BTE and find excellent agreement. We also 
demonstrate how this approach readily supports the in¬ 
clusion of detailed phonon properties, including a full 
phonon dispersion and MFP distribution. Finally, as an 
illustration of the technique, we consider a simple model 
problem using two different, commonly-used heating con¬ 
ditions. Not surprisingly, we find that the different cases 
produce different results, but we also show that the quan¬ 
tities that would be measured in an experiment are in¬ 
sensitive to the specific heating condition, at least for 
this simple, model problem. The main conclusion of this 
work is that the range of problems that can be addressed 
with diffusion equations is much broader than has been 
generally understood. 

The paper is outlined as follows. Section m presents 
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the problem under consideration and describes our the¬ 
oretical approach. Section |IIT] shows our solutions for 
the model structure, compares our solutions to those ob¬ 
tained from the phonon BTE, and applies the technique 
to bulk silicon using detailed material properties. Section 
m discusses our results and the relation to experiments. 
Finally, in Section |V] we summarize our findings. 


II. MODEL STRUCTURE AND THEORETICAL 
APPROACH 

In this work we model a simple structure comprised of 
a semi-infinite (0<a;<oo) semiconductor/insulator slab, 
driven by periodic harmonic heating at the surface (x = 
0). Actual structures in thermoreflectance experiments 
are more complicated [3, a i, and modeling such ex¬ 
periments requires considering, for example, the metal 
transducer and the finite size of the heating source in the 
y-z plane at the surface. We chose a simple structure 
to more easily illustrate and analyze the role of ballistic 
phonon effects, to demonstrate how such effects are cap¬ 
tured by diffusion equations, and to compare to rigorous 
numerical solutions of the phonon BTE [^, . 

With time-domain thermoreflectance (TDTR) the 
heating at the surface is driven by a train of short laser 
pulses ps or fs) that are modulated at a given fre¬ 
quency (/ ~ 1-10 MHz), while with frequency-domain 
thermoreflectance (FDTR) the heating at the surface is 
generated by a modulated continous laser stream (/ ^ 1- 
100 MHz). Here, we consider the simpler case of periodic 
harmonic heating, similar to the conditions of FDTR, al¬ 
though it is possible to also calculate the TDTR response 
since both are mathematically connected (^ . 

We begin with the McKelvey-Shockley flux method 
[ 2 ^ [ 2 ^ , which was shown to treat phonon transport from 
the ballistic to diffusive transport regime [ 2 ^: 
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where /g(a;,t, e) are the forward/backward heat fluxes, 
A(e) is the mean-free-path for backscattering, is 

the average ^-projected velocity, and e is the phonon en¬ 
ergy. The net heat current and heat density are given 
by Iq = Iq - Iq and Q = {1++ Iq)/v+. Eqns. 

@ have been derived assuming each phonon energy, or 
equivalently phonon frequency {v = e/h), is indepen¬ 
dent (i.e. scattering treated at the level of relaxation 
time approximation), and that the angle-dependent x- 
projected phonon velocity distribution is approximated 
by the angle-averaged value, v'^. Although our approach 
is written in terms of phonon energy e, if preferred, it 
is possible and equivalent to deal with phonon frequency 
by making the substitution e —> hiz. 


Under conditions of small temperature variations, the 
McKelvey-Shockley equations can be rewritten exactly 
as the hyperbolic heat equation (HHE) and the Cattaneo 
equation |23l |: 
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where T is the temperature, tq = A/(2u+) is the heat 
relaxation time, k = CyAu^/2 is the bulk thermal con¬ 
ductivity and Cv is the heat capacity. The equivalence 
between this expression for bulk thermal conductivity 
and the classic relation is shown in Appendix B of Ref. 
[ 2 ^. (Note that the temperature in these equations is 
the average of the temperature of the forward and re¬ 
verse heat fluxes.) These equations are modified versions 
of the heat equation and Fourier’s law that capture finite- 
velocity propagation ■ While these diffusion equations 
are widely believed to break down when ballistic effects 
are present, we showed that this is not the case 1^. 
Here we will solve the HHE and Cattaneo eq. to study 
the role of ballistic effects in frequency-dependent ther¬ 
mal transport. 

Since periodic harmonic oscillations are driving the 
heating at the surface, one can show that the solutions 
for temperature and heat current will have the form 
T{x,t) = r(x)e“* and lQ{x,t) = /Q(x)e®‘^*, respectively. 
Inserting these expressions in the HHE and Cattaneo eq., 
we obtain the following new equations: 

K rP-T 

Iq = -7T—^(6) 

(1 -I- iwtq) ox 


which can be viewed as modified versions of the heat 
equation and Fourier’s law, and as tq—^0 we retrieve 
these classic expressions. We can see how the heat cur¬ 
rent will be reduced by the denominator at high fre¬ 
quency compared to that predicted by Fourier’s law. We 
will show that by solving these simple diffusion equa¬ 
tions, we can obtain excellent agreement compared to 
numerical results of the phonon BTE without parameter 
adjustment, as shown in Fig. [2Kb)-(e) and Fig. [3] (solid 
lines: our approach, markers: phonon BTE). 


A. Boundary conditions 

Here we present the correct physical boundary con¬ 
ditions, equivalent to those from the phonon BTE, to 
be used with our approach based on diffusion equations 
(i.e. the HHE and Cattaneo eq.). Two types of peri¬ 
odic heating at the surface have been used to study bal¬ 
listic phonon effects. Case I: the surface is in contact 
with an ideal reservoir with an oscillating temperature 
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T(t) = ATe“‘, where w is the angular frequency [il,®. 
Case II: there is an oscillating heat current at the sur¬ 
face Iqit) = 0. In both cases, another boundary 

condition imposes that the excess temperature variations 
decay to zero as x—>-oo. 

We previously showed that implementing the correct 
physical boundary conditions, that is the boundary con¬ 
ditions imposed on the directed heat fluxes Iq , is key to 
capturing ballistic effects. We find the correct physical 
boundary conditions at a: = 0 are 


T{0+,t)- 


X ^ 

2(1-1- iojtq) dx 


= ATe^' 


0+ 


lQiO+,t) = I^Qe 
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(7) 
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for the case of T-controlled (case I) and Ig-controlled 
(case II) heating at the surface, respectively. Appendix 1X1 
shows how these boundary conditions are obtained. Note 
that a traditional approach would not have the second 
term on the left-hand side of Eq. 0, which is responsible 
for capturing temperature jumps at the surface. In what 
follows, “traditional approach” refers to solving the heat 
equation and Fourier’s law using the classical boundary 
conditions, i.e. r(0+,t) = or lQ{0+,t) = 


III. RESULTS 


A. Analytical solutions 


Solving the HHE and Cattaneo eq., Eqns. ©-([B]) with 
the appropriate boundary conditions, we obtain analyti¬ 
cal solutions for the temperature and heat current distri¬ 
butions. For the case of a temperature-controlled surface 
(r(t) = ATe*"^*, given by Eq. (O), we find 


T{x,t) = AT 
lQ{x,t) = AT 


1 -I- IUJTq 


(1 -I- ioJTQ) + Xk/2 
nk 


(1 -I- iujtq) + Xk/2 
where fc(a;) is expressed as 
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k{w) = ^^(l + iwrQ). (11) 

For the case of a heat current-controlled surface {Iq (t) = 
given by Eq. ([8])), we find 


Tix,t)=I°Q 

lQ{x,t) = I^e 
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The wavenumber, fc, describes the spatial distribution 
of T{x,t) and lQ{x,t). It has real and imaginary parts 
that control the decaying and oscillating components of 
the solutions, respectively, which are plotted versus fre¬ 
quency in Fig. [TJa). At low frequency (w I/tq) 




FIG. 1: (a) Real and imaginary parts of k wavevector versus 
frequency lotq. Traditional approach corresponds to solving 
the heat equation, (b)-(c) Magnitude and real parts of the 
normalized temperature profile T{x) versus normalized posi¬ 
tion x/Lp for u>tq=1Q~^ (b) and 10^ (c), where Lp=Ke{k) 
is the penetration depth. The adopted parameters are taken 
from Ref. and discussed in the caption of Fig. [21 


Re[fc] = Im[fc] = \Jio/Xvt, which can be rewritten as 

y/ ttCvZ/k the well-known classical expression for pen¬ 
etration depth At high frequency (w ^ ^/tq) 

Re[fc] = 1/A and Im[fc] = w/uj/, indicating that phonons 
on average cannot decay on a length scale shorter than 
the MFP, and that phonons travel at their ballistic ve¬ 
locity v//. 


The temperature profiles for low and high frequency 
are shown in Fig. |ljb)-(c). With the traditional ap¬ 
proach, T{x) looks like Fig. [Ijb) at all frequencies, while 
with the HHE the shape changes at higher frequency 
(Fig. HKc)). At low frequency transport is diffusive, and 
at higher frequencies transport becomes quasi-ballistic 
(purely ballistic when A —>■ oo). We find that the key 
quantity controlling the transition from diffusive to quasi- 
ballistic transport is the time tq relative to 1/w, as pre¬ 
viously highlighted by Yang and Dames [l^. It is im¬ 
portant to keep in mind that tq = X/{2v//) varies with 
dimensionality. In the case of isotropic dispersion and 
scattering time, we have A = (4/3)ugr in 3D, {TT/2)vgT 
in 2D and 2vgT in ID, as well as = Vg/2 in 3D, (2/7r)wg 
in 2D and Vg in ID, where Vg is the group velocity and r 
is the phonon scattering time. This gives tq = (4/3)t in 
3D, (7r^/8)r in 2D and t in ID. See Appendix iBl for an 
alternative formulation of tq. 
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FIG. 2: Thermal response with the temperature-controlled 
condition at the surface, (a) Normalized temperature pro¬ 
file {T(x) — To)/AT versus normalized position x/X for 
LOTQ=10~^, 10^, where To is the background temperature. 
Thick solid lines are solutions to the HHE with the physi¬ 
cally correct boundary conditions given by Eq. ((7|). T^ in¬ 
dicate the temperature of the forward/backward phonon dis¬ 
tributions, with T = (T'*' -|-r“)/2. The traditional approach 
corresponds to solving the heat equation. The surface tem¬ 
perature (b), surface heat current (c), penetration 

depth Lp (d) and phase diffenrence between and 
are plotted versus frequency lotq. Markers are results of the 
phonon lattice BTE (LBTE, taken from M)- We adopted 
the parameters in [2^ for this ID problem (see Appendix [B|: 
A=2x41 nm, u^=6733 m/s and (7^ = 1.66x10® Jm“®K~^. 


B. Comparison to the phonon Boltzmann 
transport equation 


1. Case of temperature-controlled surface 

Using Eq. Fig. UKa) shows the normalized 

temperature profile {T{x) — To)/AT (Tq is the back¬ 
ground temperature) versus normalized position (x/X) 
for ojtq=10~‘^, 10^. Solutions of our approach, using the 
HHE with the correct physical boundary conditions (Eq. 
0), are shown as thick solid lines. Solutions to the tra¬ 
ditional approach, using the heat equation (HE) with the 
traditional boundary condition T^®(0'^, t) = are 

shown as dashed lines. 


At the lower frequency both our approach and the 
traditional approach yield similar temperature profiles, 
while at higher frequency we observe significant devia¬ 
tions near the surface. At the higher wrg, and shorter 
time scales, phonons do not have sufficient time to scat¬ 
ter enough to achieve near local equilibrium. Thus, 
the phonons travel quasi-ballistically which leads to a 
nonequilibrium phonon distribution. 

To visualize the out-of-equilibrium phonon population, 
we plot the temperature profiles of the forward and back¬ 
ward moving phonons T^{x) = T{x) ± R^^'^Iq{x)/2, 
where = 1 = Cyv/f/2 is the ballistic ther¬ 
mal conductance [^, . It is important to note that 

our derivation of the HHE and Cattaneo eq. starting 
from the McKelvey-Shockley equations does not assume 
local thermal equilibrium [^. Allowing both halves of 
the phonon population to be different and to have sepa¬ 
rate temperatures is key to capturing ballistic trans¬ 
port effects. The temperature that appears in Eq. m is 
simply the average of both forward and reverse tempera- 
tures, T = (r+ -b T-)/2 When T+ and T" are 

close, the phonons are near equilibrium and transport is 
diffusive. When there is a large splitting between T+ and 
T“, the phonons are out of equilibrium and transport is 
quasi-ballistic, as seen for ujtq=10^. 

The forward-moving phonons injected at the surface 
are in equilibrium with the contact, T+(0+, t) = 
but T“(0“'') depends on how many injected phonons have 
time to scatter and return to the surface as backward- 
moving phonons. As frequency increases, this probability 
decreases, along with T“(0+). This explains the tem¬ 
perature jump observed at a:=0+, since temperature is 
the average of both streams T = (T+ -b T~)/2. This 
is mathematically equivalent to the case of an interface 
resistance equal to half the ballistic thermal resistance 
R^^^^Iq{0^)/2, although we consider ideal reflectionless 
contacts. 

In Fig. [5] (b)-(e) we present the surface temperature 
j^surf, gu^j-face heat current penetration depth Lp 

and phase difference between and , versus fre¬ 
quency UJTQ. Thick solid lines are solutions to our ap¬ 
proach and markers are numerical results of the phonon 
lattice BTE (LBTE) (Ref. [13|)- Excellent agreement is 
observed. The adopted parameters are those from Ref. 


Tradionally, as ui increases the penetration depth de¬ 
creases which drives a higher heat current (from Fourier’s 
law). We see that at higher frequency Lp^X and 
where is the ballistic heat current (the 
largest possible heat current). The phase changes from 
45° to 0° as uj increases, since only the forward-moving 
phonons contribute to the temperature and the heat cur¬ 
rent as transport becomes more ballistic. Thus 
and /g“''^ both respond instantaneously to the energy in¬ 
jected from the contact, i.e. the heat is carried away 
from the surface as efficiently as possible. The traditional 
approach breaks down at high frequency when ballistic 
effects become important, however the HHE and the Cat- 
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taneo eq. are shown to extend the traditional approach 
to much higher frequencies. 

We note that our solutions, given by Eqns. (l^- dKll) 
are identical to those reported by Yang and Dames 
for a ID problem. In fact by adding their flux equations 
derived from the BTE (Eqns. (A6)-(A7)) we directly ob¬ 
tain the HHE, thus indicating that a solution of the ID 
BTE is equivalent to that of the HHE. For a 3D material, 
where the angle-dependence of phonon transport must be 
considered, both approaches show that the ID equations 
remain valid with a rescaling of the input parameters (as 
we discussed above), which leads to some small numeri¬ 
cal differences. Our results are also identical to those re¬ 
ported by Regner et al. [l5l| if we multiply our thermal re¬ 
laxation time by 3/4, which is equivalent to t = (3/4)rQ 
the phonon scattering time. These differences originate 
in how the angle-dependence of phonon transport is ap¬ 
proximated; we replace the x-projected phonon velocity 
distribution with its angle-averaged value, v'^, while in 
Refs. [isl . solutions come from taking the two low¬ 
est order moments of the BTE. Note that an isotropic 
medium assumption is not required, and that the input 
parameters u,/ and A can be extracted for any phonon 
dispersion and scattering time (see Eqns. (IB1I) - (IB2I) '). 


2. Case of heat eurrent-controlled surface 

Using Eq. (EH), Fig. IHa) shows the temperature pro¬ 
file T{x) versus normalized position {x/X) for two fre¬ 
quencies. Solutions of our approach are shown as thick 
solid lines. Markers are numerical results of the LBTE 
(Ref. i). Excellent agreement is observed. The adopted 
parameters are those from Ref. @ , but we found that we 
had to multiply tq in Eq. ([3]) by two, which is probably 
due to a different definition of tq. A significant split¬ 
ting in T+ and T~ is observed, a signature of nonequi¬ 
librium phonons arising from ballistic effects, which be¬ 
comes more pronounced at higher frequency. 

Fig. (Hb) presents the surface temperature versus 
WTQ. Contrary to the temperature-controlled case, the 
heat current-controlled case gives a temperature that is 
larger than that expected from the traditional approach. 
The heat current can be written as Iq = — T~), 

thus the magnitude of T+ —T“ at the surface is constant. 
At higher frequencies the contribution to T~ decreases, 
since injected phonons are less likely to scatter on short 
time scales and return to the surface, and only the excess 
forward-moving phonons carry the heat current. This re¬ 
quires that /q(0“''), the magnitude of which is a constant, 
approaches Ar’’‘^*^dT+(0+) as a; where (5T+ is 

the excess temperature variation around the background 
temperature. 


Semiconductor/insulator — 



0 1 2 3 4 5 



(DTq 

FIG. 3: Thermal response with the heat current-controlled 
condition at the surface, (a) Temperature profile T{x) ver¬ 
sus normalized position x/X for /=a;/27r=7.2x 10® Hz and 
7.6x10^° Hz. (b) Surface temperature versus frequency 
uJTQ. Thick solid lines are solutions to the HHE with the phys¬ 
ically correct boundary conditions given by Eq. ®. Markers 
are results of the phonon LBTE (taken from i). r± indi¬ 
cate the temperature of the forward/backward phonon distri¬ 
butions, with T = (T'*' -I- T~)/2. The traditional approach 
corresponds to solving the heat equation. We adopted the pa¬ 
rameters in for this ID problem using effective 3D parame¬ 
ters (see Appendix IbI) : A=(4/3)q:x 40 nm, v/f=6733/{2a) m/s 
and C'v = 1.66x 10® Jm“®K~^. We find that 0 =%/^reproduces 
the numerical results of the BTE in [^. This is equivalent to 
multiplying tq in Eq. m by two, which is likely due to a 
difference in definition of tq. 


3. Apparent thermal conductivity 

It is common to define an apparent thermal conductiv¬ 
ity, Kapp, that captures the effect of non-diffusive phonon 
transport through a reduction in bulk thermal conduc¬ 
tivity as ballistic effects become prominent. Kapp is de¬ 
fined as the thermal conductivity that is extracted by 
assuming heat transport can be described by traditional 
Fourier’s law and heat eq uation. Using the definition 
Kapp = \Iq/ i—9T/dx)\ with our solutions for T{x,t) 

and lQ{x,t) we obtain 

^app — ^bulk/1 4” 

where the x and t dependences cancel out. This expres¬ 
sion is insensitive to the choice of boundary type (i.e. 
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temperature-controlled or heat current-controlled). We 
find Kapp —)> Kbuik when ojtq « 1, as expected, and 
i^app —t Kbuik /( wTQ ) when cotq >> 1. This simple equa¬ 
tion for Kapp is identical to that reported by Yang and 
Dames . although both expressions appear different. 
In Ref. [20j a frequency of interest is defined as when 
Kapp = Kbuik/2, which is numerically determined to be 
cotq = 1.73. A straightforward evaluation of Eq. (fHl) 
shows this condition corresponds to ujtq = ~ 1.73. 

Using the same definition for Kapp, our results are 
consistent with those reported by Regner et al. [H, 
given that our solutions for T[x,t) and lQ{x,t) in 3D 
are identical if we replace tq —>• r (as discussed above). 
The suppression function in this case is given hy S = 
Kapp/^buik = l/\/l -b (ujtq)^. This result suggests that 
it may be more convenient, in the case of this par¬ 
ticular model problem, to integrate over heat relax¬ 
ation time to obtain the apparent thermal conductiv¬ 
ity as opposed to the mean-free-path, i.e. Kapp(w) = 
/“ S{uj,tq) Kbuik(TQ) drg (a point highlighted by Yang 
and Dames [1^). In general, mean-free-path may be the 
more convenient integration quantity, for example when 
including the effect of laser spot size. 


C. Full phonon dispersion and mean-free-path 
distribntion: case of bulk silicon 

The results presented up to this point have been within 
the gray approximation, considering only a single phonon 
velocity and MFP. In realistic materials, the full phonon 
dispersion and energy-dependent MFP distribution must 
be treated. The most straight-forward extension is to 
consider the phonon energy (frequency) channels as in¬ 
dependent. When deriving the McKelvey-Shockley flux 
equations, with scattering treated at the level of the re¬ 
laxation time approximation, the energy channels de¬ 
couple. This is clearly an approximation and concerns 
have been raised [1^, but comparisons to full solutions 
of the phonon BTE in the steady-state [13 and transient 
[13 cases show reasonable agreement. In this section we 
demonstrate how our approach can treat realistic materi¬ 
als by using analytical solutions at each energy and then 
performing the appropriate integration over all energy 
channels. 

As a case study, we consider bulk silicon with the 
full phonon dispersion extracted from first principles cal¬ 
culations, and including boundary, defect and phonon- 
phonon Umklapp scatterings treated with phenomeno¬ 
logical models calibrated to experimental data. This 
model provides good agreement with both the experi¬ 
mental phonon energies and the measured temperature 
dependence of the thermal conductivity. Details can be 
found in Refs. Using the detailed material prop¬ 

erties of Si we extract v^{€) and A(e) (see Appendix IbI) . 
which are used to evaluate fc(e), K(e) and TQ(e) appearing 
in our solutions of T (Eq. @ and Eq. (IT^ l and Iq (Eq. 
(fTUl) and Eq. (fT^l. 
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FIG. 4: Thermal response of bulk silicon using full phonon 
dispersion and mean-free-path distribution, (a) Normalized 
temperature profile {T{x) — To)/AT versus normalized posi¬ 
tion x/X for /=a;/27r=10i^Hz, 10®Hz and 10®Hz, where To is 
the background temperature. A temperature-controlled con¬ 
dition is used, (b)-(c) Surface temperature and surface 
heat current versus frequency ujtq. Thick solid lines 

are solutions to the HHE and Cattaneo equation. Dashed 
lines correspond to the gray approximation using A=151 nm, 
u,J = 1255 m/s and (7^ = 1.63x10® Jm“®K“i. Thin solid lines 
are solutions to the heat equation and Fourier’s law. We used 
the same full dispersion and mean-free-path distribution of 
bulk silicon as Refs. [^. [^. 


We compute T{x,t,e) and lQ{x,t,e) at each energy, 
and obtain the total temperature and heat current using 


pOC 

T{x,t)= T{x,t,e)Cv{e)de/Cv, (15) 

Jo 

poc 

lQix^t)= lQ{x,t,e)de, (16) 

Jo 

where Cv{f) = eU(e) [duto/dT] is the energy-dependent 
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heat capacity, D{e) is the phonon density of states, no 
is the equilibrium Bose-Einstein distribution evaluated 
at the reference temperature Tq, and Cv = /g C'v(e)de 
is the total heat capacity. Since each phonon energy is 
assumed to be independent, our analytical solutions can 
be evaluated at each energy with the final results ob¬ 
tained from the appropriate energy integration provided 
by Eqns. (|T5]) - (|TO)) . 

Fig. mja) shows the temperature profile T{x) versus 
normalized position (a;/A) in bulk Si for / = a;/27r=10^, 
10®, 10® Hz. Full solutions to our approach (solid lines) 
are compared to solutions evaluated at a single energy us¬ 
ing the average phonon properties, i.e. the gray approxi¬ 
mation (dashed lines). Small differences between the full 
and gray results are observed at the surface; larger dif¬ 
ferences occur inside the material. This indicates that 
the gray approximation within this approach can lead to 
large errors. Another observation, most clearly seen at 
higher frequency, is that the full solutions do not yield 
exponential temperature profiles. 

Fig. |4Kb)-(c) presents the surface temperature 
and heat current versus frequency, showing both 

cases of temperature- and heat current-controlled condi¬ 
tions. We compare full solutions (thick solid lines) to the 
gray approximation (dashed lines) and the traditional ap¬ 
proach (thin solid lines). When considering surface prop¬ 
erties, the full and gray solutions are reasonably close 
and within roughly a factor of two. The surface tem¬ 
perature and heat current are sensitive to the choice of 
boundary type (T- versus Ig-controlled), and the differ¬ 
ences increase as the frequency decreases. This is most 
prominent at lower frequencies where transport is diffu¬ 
sive, and is, in fact, well-known from classical thermal 
physics. We also note that although the solutions appear 
different the thermal conductivity and thermal resistance 
are invariant to the choice of either T- and /g-controlled 
cases. Later we highlight how ballistic effects do lead to 
differences in the solutions for both cases. 


IV. DISCUSSION 

Fig. |4Kb)-(c) shows that deviations from the tradi¬ 
tional approach appear at different frequencies depending 
on the adopted heating case. We define the diffusive-to- 
ballistic transition frequency as the frequency at which 
the full solution differs from the traditional solution by 
5%. In Fig. SKb) this transition occurs at 3.6x 10^ Hz and 
3.2x10® Hz for the T- and /g-controlled cases, respec¬ 
tively, a factor of roughly 10®. Ballistic effects can come 
in through either the governing equations (i.e. Eqns. 

(g])) and/or the boundary conditions (i.e. Eqns. Q-®). 
The former cannot explain the difference in transition fre¬ 
quency, since we solve the HHE and Cattaneo eq. in both 
cases. For the /g-controlled case the boundary condition 
(Eq. ®) is the same in the traditional limit (i.e. with 
no ballistic effects), however for the T-controlled case the 
boundary condition (Eq. ([7])) differs from the traditional 


limit due to the term with dT/dx. This extra term, re¬ 
sponsible for the temperature jump at the surface, drives 
the full solution away from the traditional solution at a 
lower frequency compared to the /g-controlled case. In 
Fig. SKc) the transition frequency is 1.4x10^ Hz for the 
T-controlled case, while no deviation from the traditional 
solution is observed for the /g-controlled case (boundary 
condition for both the full and traditional solutions are 
the same). 

When analyzing experimental data, several factors 
should be considered. For example, it is not clear, which 
case (T- or /g-controlled) should be used at the surface 
of the semiconductor. We note, though, that phase lag 
is a key quantity in the analysis of FDTR and TDTR 
experiments, and that our results indicate phase is insen¬ 
sitive to the adopted heating type. One-dimensional so¬ 
lutions are also probably not adequate, and it is not clear 
if and how the metal itself and the metal-semiconductor 
junction should be treated [l^. The diffusive-to-ballistic 
transition frequencies reported here are quite high, which 
could be a result of the simplified model we used to de- 
montrate the technique. 


V. SUMMARY 

We have shown that when the correct physical bound¬ 
ary conditions are used, the HHE and Cattaneo equation 
(i.e. diffusion equations) can be used to model ballistic 
effects in frequency-dependent transient thermal trans¬ 
port. Our analytical solutions, derived for the case of 
temperature- and heat current-controlled heating at the 
surface, are found to reproduce rigorous solutions of the 
phonon BTE with high accuracy. Numerical solutions 
of the BTE were available for only ID transport prob¬ 
lems within the gray approximation, so our approach re¬ 
mains to be tested in cases where a treatment of angle- 
dependent phonon transport is required and scattering is 
handled beyond the relaxation time approximation. 

By calculating the thermal transport response of bulk 
silicon, we demonstrated how the approach can easily 
handle a full phonon dispersion and energy-dependent 
MFP distribution when the energy channels are treated 
as independent. Testing the gray approximation, we 
found that it performs resonably well at the semiconduc¬ 
tor/insulator surface, but it fails to accurately describe 
the temperature profile inside the material. 

An advantage of this simple approach based on solv¬ 
ing diffusion equations is its physical transparency. For 
example, we discussed how results can be explained in 
terms of the directed temperatures, T^, (the tempera¬ 
ture of each half of the phonon distribution) and how a 
large splitting of T'*' and T~ results from the nonequilib¬ 
rium nature of ballistic transport. The main conclusion 
of this work is that diffusion equations have the poten¬ 
tial for accurately treating ballistic to diffusive thermal 
transport. 
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Appendix A: Boundary conditions at the surface 


Here we derive the correct physical boundary condi¬ 
tions, that come from the phonon BTE, to be used when 
solving the HHE and the Cattaneo equation. Below we 
consider two cases: when heating at the surface is driven 
by a controlled temperature or a controlled heat current. 
Our starting point are the directed heat fluxes Iq used in 
the McKelvey-Shockley equations, and how to rewrite the 
boundary conditions for Iq into boundary conditions for 
T and Iq. Instead of dealing with Iq we can work with 
the directed phonon temperatures = 6Iq / -|- Tq 

[13, nil (i.e. temperature of each half of the phonon dis¬ 
tribution), where SIq is the variation in directed heat 
flux around the background equilibrium flux and Tq is 
the background temperature. 

We begin by writing two general expressions that relate 
to T and Iq [H, [111 : 


T{x,t)= [r+(x,t)+T-(x,f)] /2, 

(Al) 

lQ{x,t) = K'^^^^[T+{x,t)-T-{x,t)]. 

(A2) 

By eliminating T”, we find 


T oball 

r+(o+,t) = r(o+,t) -h ^ . 

(A3) 

Using = Cyvtl2 [H [H with the Cat- 


taneo equation, we obtain 


r(o+,t)- 


A ^ 

2{l + iojTQ) dx 


= r+(0+,t). (A4) 


0+ 


For the temperature-controlled case, where the surface of 
the semiconductor/insulator is joined to a thermal reser¬ 
voir via an ideal reflectionless contact, the temperature 
of the injected phonons is equal to the temperature of the 
reservoir. This gives the boundary condition that is Eq. 
©• For the heat current-controlled case, the HHE and 
Cattaneo equation can be solved using Iq (0+, t) = IqC^^*' 
(Eq. ([U), which states that |r+(0+, f) — r“(0+, t)| must 
be a constant. 


Appendix B: Definition of tq 


As shown in the main text and in Ref. [Ill the thermal 
relaxation time is written as tq = \/{2v^). The mean- 
free-path for backscattering and average x-projected ve¬ 
locity are defined as [H, [l^ 


A(e) = 2 


(ug(fc)r(fc)) 

{\vAm ’ 


vtie) = (bx(fc)l). 


(Bl) 

(B2) 


where r is the phonon scattering time, {X) = 

Xlfc Ai(A:)i5(e—e(fc))/Xlfeis the phonon en¬ 
ergy dispersion relation (related to the phonon frequency 
dispersion through v{k) = e{k)/h), fc is a vector in recip¬ 
rocal space, and the summation over k is restricted to 
the Brillouin zone. Using the above relations for A and 
and inserting them into our expression for tq we find 


'rQ(e) 


{vl{k)T{k)) 

{\vT{kW 


(B3) 


Thus tq is related to t through Eq. (IB3D which depends 
on the phonon dispersion. 
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